The SSU rRNA gene is one of the most widely utilized loci for phylogenetic inference among eukaryotic organisms. Although they have an average length of 1800 to 1900 bp, several unusually large 18S rDNA sequences have been reported. After examining GenBank sequences and 180 new 18S rRNA sequences from several metazoan groups, we report many other extraordinary sequences ranging between ca. 1350 bp (in symphylan myriapods) to ca. 3300 bp (in some strepsipteran insects). Myriapods are particularly interesting, having independently evolved extraordinary sequences in the four classes (Chilopoda, Diplopoda, Symphyla, and Pauropoda). An insertion event of ca. 300 bp has been detected in all but the most basal family of geophilomorphan centipedes. Other major insertions are also found in other arthropod groups, in onychophorans, molluscs, chaetognaths, echinoderms, and parasitic platyhelminths. The use of information derived from secondary structure predictions combined with a new method to analyze DNA sequence data without multiple sequence alignments is proposed as a solution for analyzing sequence data that possess alternatively conservative and variable regions, such as ribosomal genes.
INTRODUCTION
The small-subunit ribosomal RNA gene is one of the most widely utilized loci in phylogenetic inference among eukaryotic organisms (e.g., van de Peer and De Wachter, 1997) especially for the examination of phylogenetic relationships among metazoans. Ribosomal genes play a fundamental role in the synthesis of proteins in eukaryotic and prokaryotic cells. The SSU rRNA locus (in particular the 18S rRNA gene) has been widely used to infer phylogenetic relationships for reasons that have been elaborated elsewhere (e.g., Sogin, 1991; Adoutte and Philippe, 1993). The main rationale is (1) that this gene is long enough to provide phylogenetic information, (2) it is among the slowest evolving sequences found in all living organisms, and (3) different regions of the molecule evolve at different rates. For these reasons the 18S rRNA gene allows the inference of phylogenetic history across a broad taxonomic range. Moreover, the presence of many copies per genome and their homogenization through concerted evolution (Dover, 1982; Hillis and Dixon, 1991) greatly reduce intraspecific variation and facilitate DNA amplification via PCR.
The SSU rRNA gene of most organisms is between 1800 and 1900 bp in total length. However, exceptional cases of long SSU rRNA genes are not uncommon (see Crease and Colbourne, 1998 for some examples), with certain regions being particularly prone to sequence variability (both in nucleotide substitutions and in sequence length).
In some cases, the SSU rRNA gene contains type I introns that are spliced out of the mature rRNA molecule (e.g., De Wachter et al., 1992; Wilcox et al., 1992; Bhattacharya et al., 1994). However, there are many long SSU rRNA genes that do not contain introns (e.g., Hinkle et al., 1994). Some of these long SSU rRNA genes have been reported for metazoans: Acyrthosiphon pisum 2469 bp (Kwon et al., 1991); Xenos vesparum 3316 bp (Chalwatzis et al., 1995); Daphnia pulex 2293 bp (Crease and Colbourne, 1998); Euperipatoides leuckarti 2206 bp (Aguinaldo et al., 1997); Echinococcus granulosus 2394 bp (Picón et al., 1996). To our knowledge, type I introns have not been found in any metazoan taxon.
Other unusual phenomena have been reported for some organisms. The presence of more than one type of SSU rRNA gene has been found in the protozoan Plasmodium (Gunderson et al., 1987; Waters et al., 1989; Qari et al., 1994). Two types of 18S rRNA genes have also been reported in a group of metazoans, the freshwater and terrestrial planarians of the family Dugesiidae (Carranza et al., 1996, 1998a, 1998b).
Because alternate conserved and nonconserved regions are present in the SSU rRNA, sequence alignments are challenging. Some authors have based alignments on information derived from secondary structure models (see Kjer, 1995). However, different models (i.e., Gutell et al., 1985; Hendriks et al., 1988a, 1988b; Neefs and De Wachter, 1990; Van de Peer et al., 1998) can lead to different phylogenetic results (Winnepenninckx and Backeljau, 1996). Moreover, the variable regions cannot be aligned reliably and thus used as phylogenetic information, even if secondary structure predicts that a certain string of nucleotides is homologous (e.g., a variable loop, situated between a conserved stem).
The existence of different models can be attributed to problems in inferring secondary (and tertiary and quaternary) structures via direct observation. This was pointed out by Gutell and collaborators, who acknowledged that “any rigorous search for a secondary structure model for 16-S rRNA would necessitate use of the comparative method” (Gutell et al., 1985: 156). These authors also mentioned that the direct approach of X-ray crystallography remains only a remote possibility, because of the difficult procedure for preparing high-quality crystals of ribosomes or their subunits, together with the added problem that tertiary structure may be directed and/or stabilized by quaternary interactions (Gutell et al., 1985).
Where direct approaches have failed to infer the secondary structure of the SSU rRNA gene, some indirect approaches have succeeded. The comparative method has been successful to the point of building a database of hundreds of secondary structures of SSU rRNA. This database, the SSU ribosomal subunit RNA database ( http://rrna.uia.ac.be/rrna/ssu [Van de Peer et al., 1998]) is maintained and updated constantly, and incorporates probably the largest comparative molecular data set.
In the present study, we report some extraordinary examples showing enormous variation in length among 18S rRNA sequences of different metazoan taxa. We also comment on the regions of the molecule that concentrate the largest variation. Finally we explore a novel method that facilitates the use of variable regions of the molecule, which cannot be accommodated into standard phylogenetic analyses using sequence alignments.
Materials and Methods
Genomic DNA samples were obtained from fresh, frozen, or ethanol-preserved tissues in a solution of Guanidinium thiocyanate homogenization buffer following a modified protocol for RNA extraction (Chirgwin et al., 1979). The 18S rDNA locus was PCR-amplified in two or three overlapping fragments of about 950, 900, and 850 bp each, using primer pairs 1F–5R, 3F–18Sbi, and 5F–9R, respectively. Primers used in amplification and sequencing were described in Giribet et al. (1996, 1999). Amplification was carried out in a 50 μL volume reaction, with 1.25 units of AmpliTaq DNA Polymerase (Perkin Elmer), 200 μM of dNTPs, and 1 μM of each primer. The PCR program consisted of a initial denaturing step at 94°C for 60 seconds, 35 amplification cycles (94°C for 15 sec, 49°C for 15 sec, 72°C for 15 sec), and a final step at 72°C for 6 minutes in a GeneAmp PCR System 9700 (Perkin Elmer).
PCR samples were purified with the GENECLEAN III kit (BIO 101 Inc.) and directly sequenced using an automated ABI Prism 377 DNA sequencer. Cycle-sequencing with AmpliTaq DNA Polymerase, FS (Perkin-Elmer) using dye-labeled terminators (ABI PRISM BigDye Terminator Cycle Sequencing Ready Reaction Kit) was performed in a GeneAmp PCR System 9700 (Perkin Elmer). Amplification was carried out in a 10 μL volume reaction: 4 μL of Terminator Ready Reaction Mix, 10–30 ng/mL of PCR product, 5 pmoles of primer and dH20 to 10 μL. The cycle-sequencing program consisted of a step at 94°C for 3 minutes, 25 sequencing cycles (94°C for 10 sec, 50°C for 5 sec, 60°C for 4 min) and a rapid thermal ramp to 4°C and hold. The BigDye-labeled PCR products were isopropanol-precipitated following manufacturer protocol.
Some PCR products that could not be sequenced directly were purified and ligated into pUC 18 Sma I/BAP dephosphorylated vector using the SureClone Ligation Kit (Pharmacia P-L Biochemicals) as described in Giribet et al. (1996). Sequencing was then performed by the dideoxy termination method (Sanger et al., 1977) using T7 DNA polymerase (Sequencing Kit from Pharmacia Biotech).
All sequences have been deposited in GenBank (see taxonomy, sequence length, and accession codes in table 1). The new sequences have been compared to other published sequences available from GenBank. The terminology used for the secondary structure topology follows the nomenclature of Van de Peer et al. (1998).
RESULTS AND DISCUSSION
Extraordinary 18S rRNA Sequences of Metazoans
After studying the 18S rRNA gene of about 400 metazoan taxa (180 of these sequences collected by the authors; table 1), we have observed that several animal taxa possess large insertions at different regions of the molecule. This variation is shown in figure 1 and table 1. For example, within the phylum Mollusca, the cephalopods Loligo pealei (squid), Sepia elegans (cuttlefish), and Nautilus scrobiculatus present large insertions in regions V2, V4, V7, and V9. Other groups of molluscs such as the anomalodesmatan bivalves present large insertions in region V7, and some Archaeogastropoda have insertions in regions V2 and V4. But insertions are not restricted to the phylum Mollusca. Sea cucumbers (Echinodermata, Holothuroidea) and arrow-worms (Chaetognatha) present insertions in region V4; some leeches (Annelida, Hirudinea) present insertions in region V7; some parasitic planarians (Platyhelminthes) present insertions both in regions V4 and V7; and velvet worms (Onychophora) present insertions in regions V2, 11, E23–7, V7, and V9. Within the arthropods, there are many extraordinary cases among insects (see the reported strepsipteran sequences by Chalwatzis et al., 1995; 1996; Whiting et al., 1997) and certain crustacean groups. But the most bizarre case within arthropods (and perhaps for the entire Metazoa) are myriapods.
The 18S rRNA Gene of the Myriapods
Myriapods comprise four groups of terrestrial arthropods. Centipedes (class Chilopoda) and millipedes (class Diplopoda) are the two principal classes of myriapods. The two other classes of myriapods lack common names: Symphyla and Pauropoda. Prior to the analysis of Edgecombe et al. (1999), the only complete 18S rRNA sequence data for myriapods available at GenBank were eight centipedes and two millipedes (Giribet et al., 1996, 1999; Giribet and Ribera, 1998). Two centipede species of the order Geophilomorpha (Clinopodes poseidonis and Pseudohimantarium mediterraneum) present an insertion of about 300 bp at region V7, whereas all the other available sequences are fairly conserved in terms of primary sequence. However, a wider ongoing study on the 18S rRNA gene of myriapods suggests that they constitute one of the most interesting cases of 18S rRNA variation in any metazoan group.
Within the centipedes, the members of the order Geophilomorpha (15 species studied belonging to 9 families), excluding two species of the most basal family Mecistocephalidae (Edgecombe et al., 1999), exhibit insertions of about 300 bp in the region V7. This is an unusual example that shows exactly when the insertion occurred during the phylogenetic process, and illustrates the putative information of such insertions (fig. 2).
Within the millipedes, members of the family Polyzonidae display sequences longer than 2700 bp. Data for one species of Pauropoda show that the 18S rRNA is ca. 2200 bp, with several small insertions (Giribet, 1997; Giribet and Ribera, 2000). But perhaps the most unusual case among metazoan 18S rRNA sequences is the Class Symphyla. Amplification of the 18S rRNA loci of three species belonging to two genera (two species of Scutigerella from northeastern Spain and the Canary Islands, respectively, and one species of Hanseniella from Australia) yielded a product band size of about 1350 bp. Sequencing this fragment suggests a deletion of about 500 bp in the central region of the molecule.
Although it might be conjectured that in this case a nonfunctional pseudogene has been sequenced, as occurred with the 18S rRNA locus of the platyhelminth Dugesia mediterranea (Carranza et al., 1996) and other dugesiids (Carranza et al., 1998a, 1998b), this seems improbable for several reasons. First, this sequence has been obtained from three different species and in two independent laboratories. Second, none of the highly conserved primers from the “deleted” region (forward primers 4F, 18Sa0.7, 18Sa0.79, 18Sa1.0; reverse primers 4R, 18Sb5.0, 18Sb3.9, 18Sb3.0 [Giribet et al., 1996; Whiting et al., 1997]) amplified any DNA fragment when combined with primers from other regions. If the 1350 bp fragment was a pseudogene, we would expect to amplify fragments of the original gene when using the conserved primers located within the “deleted” region. Third, phylogenetic analyses including the symphylan sequences show that symphylans are arthropods related to other myriapods (Giribet, 1997; fig. 3). Fourth, amplification of DNA from an RNA source, as described by Carranza et al. (1998a), yielded a product band size of approximately 1350 bp, as expected. These facts demonstrate that a deletion of ca. 500 bp occurred in the common ancestor of these three symphylan species.
18S rRNA Variation in Metazoans
It seems that large insertions and deletions are not as constrained as was previously thought (e.g., Crease and Colbourne, 1998). These events occur in many metazoan taxa, and are commonly in regions V2, V4, V7, and V9. Other parts of region 23 (that includes the region V4) are also variable. In a phylogenetic study of about 150 arthropod 18S rRNA sequences (Giribet and Ribera, 2000), insertions at region E23–7 were observed in Onychophora and in Pauropoda while insertions at region E23–8 were observed in the pauropod species. Other insertions observed in particular taxa occur at sites 8 (in the millipede Polyzonium), 11 (in Onychophora), 29 (in Pauropoda), and 46 (in the proturan Acerentulus traeghardi). However, we only obtained sequences from one pauropod and one proturan, hence these results cannot be generalized to other members of such groups.
Certain taxa present insertions in variable regions whereas in the remaining regions the primary sequence may be conserved. For example, certain geophilomorph centipedes exhibit a small insertion at region V2 (between 10 and 80 bp compared to other centipedes) and a large insertion (about 300 bp) at region V7, whereas the remaining positions are conserved with respect to other centipedes. Other taxa not only present insertions in the variable regions, but also in the primary sequence. This is the case in the cephalopods, which have insertions in regions V2, V4, V7, and V9, and differ considerably from other molluscs in the primary sequence of the remaining regions.
Although several metazoans present insertions in the 18S rRNA, reduction of the 18S rRNA gene appears to be a rare event in evolution. To our knowledge, there are no other reported cases of 18S rDNA sequences with a deletion of the magnitude of that observed in the symphylans (ca. 500 bp). The deletion corresponds to the central region of the molecule (approximately from region 14 to E23–9). The remaining primary sequence is fairly conserved, which we assume may be functional. However reconstruction of a global secondary structure cannot be conducted using the comparative method. Another case of sequence length reduction in the 18S rRNA locus occurs in the Dicyemid mesozoans (three species of the genus Dicyema), with a total gene sequence of about 1670 bp that presents two major deletions at the variable sites V2 and V9.
The geophilomorphan centipedes that present the insertion of about 300 bp at region V7 (13 species) also display a large insertion (about 300 bp) at the D3 expansion fragment of the large subunit rRNA locus (28S rRNA). Neither of the insertions at region V7 of the 18S rRNA or at the D3 expansion fragment of the 28S rRNA locus have been found in the putative most basal geophilomorph family Mecistocephalidae (Giribet et al., 1999; Edgecombe et al., 1999). This apparent correlation between the insertions at region V7 of the 18S rRNA locus and at the D3 expansion fragment of the 28S rRNA locus (also observed in the cephalopod Loligo) could suggest a possible interaction between these subunits in the ribosome.
18S rRNA Variable Sites and Phylogenetic Analyses
The presence of certain variable sites in the 18S rRNA molecule can hinder phylogenetic analyses at the alignment step, a fact that is promoted by several researchers who avoid the use of ribosomal genes for drawing inferences (e.g., Ayala et al., 1998). In general, researchers using ribosomal genes exclude the variable regions from their phylogenetic inference step because of uncertainties in alignments (e.g., Giribet et al., 1996, 1999). But since data removal from phylogenetic analyses is also problematic (see Gatesy et al., 1993), and automatic alignments are explicit, other researchers prefer the use of automatic alignments exploring different cost matrices (e.g., Wheeler, 1995). Regardless of which of these options is the best (philosophically, computationally, or practically), the extreme length variation of some of the SSU rRNA helices may constitute a serious problem in the phylogenetic analysis of certain data sets.
A clear example is illustrated in figure 4, where we show the variable region V7 of 17 centipede taxa (see table 2 for the taxonomy). The fragment as a whole can be considered homologous because it is located between two conserved regions that constitute a stem. The most common sequence length for the V7 region in centipedes ranges between 65 and 70 bp (character state present in the five recognized orders of centipedes). However three taxa (Scolopendra, Ethmostigmus and Alipes; family Scolopendridae) display sequences between 93 and 94 bp, and four taxa (Pseudohimantarium, Henia, Clinopodes and Zelanion; belonging to four families of Geophilomorpha) exhibit sequences between 354 and 384 bp. These two clades are well defined morphologically and could also be characterized in terms of sequence length or perhaps secondary structure topology. However, a standard phylogenetic analysis of this fragment cannot be conducted successfully when base-to-base homology is required, as is the case with multiple sequence alignments. This situation is frustrating since the sequence data clearly display historical information that cannot be used phylogenetically.
A new method to analyze DNA sequence data that does not require base-to-base correspondences was recently developed by Wheeler (1999) and is discussed in greater detail there. Briefly, the method, named “fixed character states”, optimizes DNA sequence data without employing multiple sequence alignments by treating entire homologous stretches of sequence data as characters. The set of specific sequences exhibited by the terminal taxa constitutes the character states. Thus the number of states is equal to the number of unique sequences (or homologous fragments) exhibited by the data. In the example illustrated here, there is one character (region V7) with 17 states (as many as different taxa). Other situations could arise where the number of states would be smaller than the number of taxa if two or more were to share identical sequences. The salient feature of this method is the treatment of length variation. Since the sequence variation is expressed through a series of transformations between states, indel or “gap” variation only occurs as transformations between sequences, not globally among all the sequence data. This has the effect of moderating the difficulties presented by extreme nucleotide length variation at the expense of treating strings of bases as character states instead of individual nucleotides.
A matrix of transformation costs is created to relate the states to one another. The cells of this matrix are defined as the minimum transformation cost required between each pair of states based on insertion-deletion and base substitution costs (as in the calculation of an alignment score). The next operation uses this transformation matrix to diagnose a specific phylogenetic topology by means of existing dynamic programming techniques (Sankoff and Rousseau, 1975) with the number of states greatly expanded. This method has been implemented in the computer program POY (Gladstein and Wheeler, 1997) specifying the option -fixedstates (available via anonymous ftp at the site ftp.amnh.org /pub/molecular/poy/).
To illustrate how the method works empirically, we have analyzed the 17 sequences presented in figure 4 (and table 2) using the fixed character states method. The tree obtained (fig. 5) shows some lack of resolution, but also shows certain clades highly consistent with the current morphology of the group, as well as with the molecular analyses of Giribet et al. (1999) and Edgecombe et al. (1999). Scolopendromorpha is recognized as a clade that in turn includes a monophyletic clade, the family Scolopendridae, presenting an insertion of about 25 bp. Another clade is defined by the insertion of about 300 bp that groups all geophilomorph species except the mecistocephalids (Mecistocephalus and Nodocephalus the most basal group that lacks the insertion). This is encouraging considering that this topology corresponds to the analysis of just a few bases from the variable region V7. Thus, this method facilitates use of all the information (variable and conserved) from ribosomal genes.
CONCLUSIONS
Long SSU rRNA appears to be more common than claimed by some authors (e.g., Crease and Colbourne, 1998) based on its occurrence in at least seven metazoan phyla: Platyhelminthes, Mollusca, Onychophora, Arthropoda, Chaetognatha, Echinodermata, and Mesozoa. However, large deletions appear to be rare, having so far only been found in one group of arthropods (Symphyla) and in mesozoans. Probably many more taxa display extraordinary 18S rRNA genes, but this will not be discovered until sampling within each phylum increases. Nonetheless, the existence of variable regions should not discourage the use of ribosomal genes in phylogenetic analyses, especially when secondary structure predictions are combined with novel methods of DNA sequence data analysis. In this sense, the characterization of secondary structural features by means of the comparative method, and the use of these features (homologous regions) as characters with multiple states provides a powerful approach for the analysis of such data using the fixed character states method. Maybe it is at such levels that secondary structure information can best contribute to phylogenetic analyses.
Acknowledgments
We thank Jerry Harasewych for discussion and Lorenzo Prendini and Michael Whiting for many helpful suggestions. G. G. was supported by a Lerner Gray research fellowship from the AMNH.
REFERENCES
Fig. 1.
Schematic representation of the 18S rRNA locus. The gray squares represent the variable regions V2, V4, V7, and V9 with insertions (V2: Onychophora, Geophilomorpha, Cephalopoda, Archaeogastropoda; V4: Hexapoda, Crustacea, Pauropoda, Holothuroidea, Chaetognatha, Platyhelminthes, Cephalopoda; V7: Onychophora, Hexapoda, Crustacea, Pauropoda, Chilopoda, Platyhelminthes, Hirudinea, Cephalopoda, Gastropoda; V9: Onychophora, Crustacea, Cephalopoda). The black arrowheads represent particular insertions (10: Pauropoda; 11: Onychophora; E23–7: Onychophora and Pauropoda; E23–8: Pauropoda; 29: Pauropoda; 46: Protura). The black bar represents the 500 bp deletion of the Symphyla
Fig. 2.
Phylogenetic tree of the centipedes based on the combined analysis of Edgecombe et al. (1999). The arrow indicates where the insertion of ca. 300 bp at region V7 occurred during the evolution of centipedes
Fig. 3.
Phylogenetic tree based on 18S rRNA sequence data indicating the position of two symphylans (box) with respect to other myriapods (underlined taxa) in a phylogenetic analysis of arthropods (from Giribet, 1997). The two symphylans appear related to other myriapods
Fig. 5.
Phylogenetic analysis of the data from fig. 4 using the “fixed character states” method of Wheeler (1999) implemented in the computer program POY (Gladstein and Wheeler, 1997). Commands: poy -fixedstates -noleading -norandomizeoutgroup -gap 1 -maxtrees 20 -multibuild 10 -seed−1 -slop 2 -checkslop 5. The two circles illustrate the insertions of the Geophilomorpha (ca. 300 bp), and the Scolopendridae (ca. 25 bp)
TABLE 1
List of the 180 Species of 18S rRNA Sequesnces Generated by the Authors (with GenBank accession codes and sequence length [excluding a total of 46 bp from the external primers 1F and 9R];asterisks refer to noncomplete sequences, and thus the length is not reported)
TABLE 1–(Continued)
TABLE 2
Taxonomy of the Centipede Species (Chilopoda) Represented in Figures 4 and 5
